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Abstract 

The most popular approaches for the simulation of dynamic systems in computer graphics are force based. Internal 
and external forces are accumulated from which accelerations are computed based on Newton's second law of 
motion. A time integration method is then used to update the velocities and finally the positions of the object. 
A few simulation methods (most rigid body simulators) use impulse based dynamics and directly manipulate 
velocities. In this paper we present an approach which omits the velocity layer as well and immediately works 
on the positions. The main advantage of a position based approach is its controllability. Overshooting problems 
of explicit integration schemes in force based systems can be avoided. In addition, collision constraints can be 
handled easily and penetrations can be resolved completely by projecting points to valid locations. We have used 
the approach to build a real time cloth simulator which is part of a physics software library for games. This 
application demonstrates the strengths and benefits of the method. 

Categories and Subject Descriptors (according to ACM CCS): 1.3.5 [Computer Graphics]: Computational Geometry 
and Object ModelingPhysically Based Modeling; 1.3.7 [Computer Graphics]: Three-Dimensional Graphics and 
RealismAnimation and Virtual Reality 



1. Introduction 

Research in the field of physically based animation in com- 
puter graphics is concerned with finding new methods for 
the simulation of physical phenomena such as the dynamics 
of rigid bodies, deformable objects or fluid flow. In contrast 
to computational sciences where the main focus is on accu- 
racy, the main issues here are stability, robustness and speed 
while the results should remain visually plausible. There- 
fore, existing methods from computational sciences can not 
be adopted one to one. In fact, the main justification for 
doing research on physically based simulation in computer 
graphics is to come up with specialized methods, tailored to 
the particular needs in the field. The method we present falls 
into this category. 

The traditional approach to simulating dynamic objects 
has been to work with forces. At the beginning of each time 
step, internal and external forces are accumulated. Examples 
of internal forces are elastic forces in deformable objects or 
viscosity and pressure forces in fluids. Gravity and collision 
forces are examples of external forces. Newton's second law 
of motion relates forces to accelerations via the mass. So us- 



ing the density or lumped masses of vertices, the forces are 
transformed into accelerations. Any time integration scheme 
can then be used to first compute the velocities from the ac- 
celerations and then the positions from the velocities. Some 
approaches use impulses instead of forces to control the an- 
imation. Because impulses directly change velocities, one 
level of integration can be skipped. 

In computer graphics and especially in computer games 
it is often desirable to have direct control over positions of 
objects or vertices of a mesh. The user might want to attach 
a vertex to a kinematic object or make sure the vertex always 
stays outside a colliding object. The method we propose here 
works directly on positions which makes such manipula- 
tions easy. In addition, with the position based approach it is 
possible to control the integration directly thereby avoiding 
overshooting and energy gain problems in connection with 
explicit integration. So the main features and advantages of 
position based dynamics are 

• Position based simulation gives control over explicit inte- 
gration and removes the typical instability problems. 
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Figure 1: A known deformation benchmark test, applied here to a cloth character under pressure. 



• Positions of vertices and parts of objects can directly be 
manipulated during the simulation. 

• The formulation we propose allows the handling of gen- 
eral constraints in the position based setting. 

• The explicit position based solver is easy to understand 
and implement. 

2. Related Work 

The recent state of the art report [NMK*05] gives a good 
overview of the methods used in computer graphics to simu- 
late deformable objects, e.g. mass-spring systems, the finite 
element method or finite difference approaches. Apart from 
the citation of [MHTG05], position based dynamics does not 
appear in this survey. However, parts of the position based 
approach have appeared in various papers without naming it 
explicitly and without defining a complete framework. 

Jakobsen [JakOl] built his Fysix engine on a position 
based approach. His central idea was to use a Verlet inte- 
grator and manipulate positions directly. Because velocities 
are implicitly stored by current and the previous positions, 
the velocities are implicitly updated by the position manip- 
ulation. While he focused mainly on distance constraints, 
he only gave vague hints on how more general constraints 
could be handled. In this paper we present a fully general 
approach which handles general constraints. We also focus 
on the important issue of conservation of linear and angu- 
lar momenta by position projection. We work with explicit 
velocities instead of storing previous positions which makes 
damping and friction simulation much easier. 

Desbrun [DSB99] and Provot [Pro95] use constraint pro- 
jection in mass spring systems to prevent springs from over- 
stretching. In contrast to a full position based approach, pro- 
jection is only used as a polishing process for those springs 
that are stretched too much and not as the basic simulation 
method. 

Bridson et al. use a traditional force based approach for 



cloth simulation [BFA02] and combine it with a geomet- 
ric collision resolving algorithm based on positions to make 
sure that the collision resolving impulses are kept within sta- 
ble bounds. The same holds for the kinematical collision cor- 
rection step proposed by Volino et al. [VCMT95]. 

A position based approach has been used by Clavet et 
al. [CBP05] to simulate viscoelastic fluids. Their approach 
is not fully position based because the time step appears in 
various places of their position projections. Thus, the inte- 
gration is only conditionally stable as regular explicit inte- 
gration. 

Miiller et al. [MHTG05] simulate deformable objects 
by moving points towards certain goal positions which are 
found by matching the rest state to the current state of the 
object. Their integration method is the closest to the one we 
propose here. They only treat one specialized global con- 
straint and, therefore, do not need a position solver. 

Fedor [Fed05] uses Jakobsen's approach to simulate char- 
acters in games. His method is tuned to the particular prob- 
lem of simulating human characters. He uses several skeletal 
representations and keeps them in sync via projections. 

Faure [Fau98] uses a Verlet integration scheme by modi- 
fying the positions rather than the velocities. New positions 
are computed by linearizing the constraints while we work 
with the non linear constraint functions directly. 

We define general constraints via a constraint function 
as [BW98] and [THMG04]. Instead of computing forces as 
the derivative of a constraint function energy, we directly 
solve for the equilibrium configuration and project positions. 
With our method we derive a bending term for cloth which 
is similar to the one proposed in [GHDS03] and [BMF03] 
but adopted to the point based approach. 

In Section 4 we use the position based dynamics approach 
for the simulation of cloth. Cloth simulation has been an ac- 
tive research field in computer graphics in recent years. In- 
stead of citing the key papers of the field individually we 
refer the reader to [NMK*05] for a comprehensive survey. 
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3. Position Based Simulation 

In this section we will formulate the general position based 
approach. With cloth simulation, we will give a particular 
application of the method in the subsequent and in the results 
section. We consider a three dimensional world. However, 
the approach works equally well in two dimensions. 

3.1. Algorithm Overview 

We represent a dynamic object by a set of N vertices and M 
constraints. A vertex ie [1, .. . , N] has a mass mj, a position 
x; and a velocity Vj. 

A constraint j € [1,. . . ,M] consists of 

• a cardinality nj, 

• a function Cj : M 3 "' -> K, 

• a set of indices {ii ,.../„.}, £ [1, . . .N], 

• a stiffness parameter kj 6 [0. . . 1] and 

• a type of either equality or inequality. 

Constraint j with type equality is satisfied if 
Cy(x,-, , . . . ,x; n .) = 0. If its type is inequality then it is 
satisfied if Cy(x,-, , . . . ,x, n ) > 0. The stiffness parameter kj 
defines the strength of the constraint in a range from zero to 
one. 

Based on this data and a time step At, the dynamic object 
is simulated as follows: 

(1) forall vertices i 

(2) initialize x, = nS, Vj = v?,Wj = 

(3) endfor 

(4) loop 

(5) forall vertices i do v; <— v; + Atwji ext (x,-) 

(6) dampVelocities(vi ..... v^r) 

(7) forall vertices i do p, «— x; + Arv,- 

(8) forall vertices i do generateCollisionConstraints(x; — > p,) 

(9) loop solverlterations times 

(10) projectConstraints(Ci , . . . , C m +m c „„ , Pi, • • - , Pat) 

(11) endloop 

(12) forall vertices i 

(13) v ; ^( Pl --x,0/Af 

(14) Xi <- p, 

(15) endfor 

(16) velocityUpdate(vi , . . . , vjv) 

(17) endloop 

Lines (l)-(3) just initialize the state variables. The core 
idea of position based dynamics is shown in lines (7), (9)- 
(11) and (13)-(14). In line (7), estimates p,- for new locations 
of the vertices are computed using an explicit Euler inte- 
gration step. The iterative solver (9)-(ll) manipulates these 
position estimates such that they satisfy the constraints. It 
does this by repeatedly project each constraint in a Gauss- 
Seidel type fashion (see Section 3.2). In steps (13) and (14), 
the positions of the vertices are moved to the optimized es- 
timates and the velocities are updated accordingly. This is 



in exact correspondence with a Verlet integration step and 
a modification of the current position [JakOl], because the 
Verlet method stores the velocity implicitly as the difference 
between the current and the last position. However, working 
with velocities allows for a more intuitive way of manipulat- 
ing them. 

The velocities are manipulated in line (5), (6) and (16). 
Line (5) allows to hook up external forces to the system if 
some of the forces cannot be converted to positional con- 
straints. We only use it to add gravity to the system in which 
case the line becomes v,- <— v; + Atg, where g is the gravita- 
tional acceleration. In line (6), the velocities can be damped 
if this is necessary. In Section 3.5 we show how to add global 
damping without influencing the rigid body modes of the 
object. Finally, in line (16), the velocities of colliding ver- 
tices are modified according to friction and restitution coef- 
ficients. 

The given constraints C\ , . . . , Cm are fixed throughout the 
simulation. In addition to these constraints, line (8) generates 
the M co n collision constraints which change from time step 
to time step. The projection step in line (10) considers both, 
the fixed and the collision constraints. 

The scheme is unconditionally stable. This is because the 
integration steps (13) and (14) do not extrapolate blindly 
into the future as traditional explicit schemes do but move 
the vertices to a physically valid configuration p, computed 
by the constraint solver. The only possible source for insta- 
bilities is the solver itself which uses the Newton-Raphson 
method to solve for valid positions (see Section 3.3). How- 
ever, its stability does not depend on the time step size but 
on the shape of the constraint functions. 

The integration does not fall clearly into the category 
of implicit or explicit schemes. If only one solver iteration 
is performed per time step, it looks more like an explicit 
scheme. By increasing the number of iterations, however, a 
constrained system can be made arbitrarily stiff and the al- 
gorithm behaves more like an implicit scheme. Increasing 
the number of iterations shifts the bottleneck from collision 
detection to the solver. 

3.2. The Solver 

The input to the solver are the M + M co u constraints and 
the estimates p; , . . . , p^y for the new locations of the points. 
The solver tries to modify the estimates such that they sat- 
isfy all the constraints. The resulting system of equations 
is non-linear. Even a simple distance constraint C(pi,p2) = 
|Pl — P2 1 ~d yields a non-linear equation. In addition, the 
constraints of type inequality yield inequalities. To solve 
such a general set of equations and inequalities, we use a 
Gauss-Seidel-type iteration. The original Gauss-Seidel al- 
gorithm (GS) can only handle linear system. The part we 
borrow from GS is the idea of solving each constraint inde- 
pendently one after the other. However, in contrast to GS, 
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solving a constraint is a non linear operation. We repeat- 
edly iterate through all the constraints and project the par- 
ticles to valid locations with respect to the given constraint 
alone. In contrast to a Jacobi-type iteration, modifications to 
point locations immediately get visible to the process. This 
speeds up convergence significantly because pressure waves 
can propagate through the material in a single solver step, an 
effect which is dependent on the order in which constraints 
are solved. In over-constrained situations, the process can 
lead to oscillations if the order is not kept constant. 



Ap 2 




Figure 2: Projection of the constraint C(pi,p2) = |pi — 
P2I — d. The corrections Ap, are weighted according to the 
inverse masses w,- = 1/m,-. 



3.3. Constraint Projection 

Projecting a set of points according to a constraint means 
moving the points such that they satisfy the constraint. The 
most important issue in connection with moving points di- 
rectly inside a simulation loop is the conservation of linear 
and angular momentum. Let Ap,- be the displacement of ver- 
tex i by the projection. Linear momentum is conserved if 



£m,Ap, =0, 



(1) 



which amounts to conserving the center of mass. Angular 
momentum is conserved if 



£r,-xm,Ap,- = 0, 



(2) 



where the r,- are the distances of the p,- to an arbitrary com- 
mon rotation center. If a projection violates one of these con- 
straints, it introduces so called ghost forces which act like ex- 
ternal forces dragging and rotation the object. However, only 
internal constraints need to conserve the momenta. Collision 
or attachment constraints are allowed to have global effects 
on the object. 

The method we propose for constraint projection con- 
serves both momenta for internal constraints. Again, the 
point based approach is more direct in that we can directly 
use the constraint function while force based methods derive 
forces via an energy term (see [BW98, THMG04]). Let us 
look at a constraint with cardinality n on the points pi , . . . , p„ 
with constraint function C and stiffness k. We let p be the 
concatenation [p\ , . . . .p,^] 7 ". For internal constraints, C is 
independent of rigid body modes, i.e. translation and rota- 
tion. This means that rotating or translating the points does 
not change the value of the constraint function. Therefore, 
the gradient V p C is perpendicular to rigid body modes be- 
cause it is the direction of maximal change. If the correction 
Ap is chosen to be along VC P both momenta are automati- 
cally conserved if all masses are equal (we handle different 
masses later). Given p we want to find a correction Ap such 
that C(p + Ap) = 0. This equation can be approximated by 



C(p + Ap) « C(p) + VpC(p) ■ Ap = 0. 



(3) 



Restricting Ap to be in the direction of V p C means choos- 
ing a scalar X such that 



Substituting Eq. (4) into Eq. (3), solving for X and substi- 
tuting it back into Eq. (4) yields the final formula for Ap 



A P : 



C(p) 



fVpC(p) 



(5) 



|V p C(p)p 

which is a regular Newton-Raphson step for the iterative so- 
lution of the non-linear equation given by a single constraint 
For the correction of an individual point p, we have 



Ap, = -s V p ,.C(pi,...,p„), 
where the scaling factor 

C(pi,...,p„) 



(6) 



(7) 



E^|v Pj .c(pi ) ... ) p„)|2 

is the same for all points. If the points have individual 
masses, we weight the corrections Ap,- by the inverse masses 
Wj = 1 /mi, In this case a point with infinite mass, i.e. w/ = 0, 
does not move for example as expected. Now Eq. (4) is re- 
placed by 

Ap,- = Aw,-V p .C(p) yielding 



_ C(pi,...,p„) 

Z.j w j\Vpj C fah--;Vn)\ 2 

for the scaling factor and for the final correction 
Ap,- = -5w,-V p .C(pi,...,p„). 



(8) 



(9) 



To give an example, let us consider the distance constraint 
function C(pi,p2) = |pi — P2I — d. The derivative with re- 
spect to the points are V pi C(pi ,P2) = n and V p ,C(pi ,P2) = 

n with n = , Pl P2 , . The scaling factor s is, thus, s = 
IP1-P2I b 

P^, - and the final corrections 



Api = 

Ap2 = + 



W] 



w\ +W2 

W2 
W\ +W2 



| Pl -p 2 |-rf)^ — ^- (10) 
|Pi -P2I 

|Pl-P2|- A Pl ~ P2 , (ID 



IP1-P2I 

which are the formulas proposed in [JakOl] for the projec- 
tion of distance constraints (see Figure 2). They pop up as a 
special case of the general constraint projection method. 



Ap = AV p C(p). 



(4) 
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We have not considered the type and the stiffness k of 
the constraint so far. Type handling is straight forward. If 
the type is equality we always perform a projection. If 
the type is inequality, the projection is only performed if 
C(pi , . . . , p„) < 0. There are several ways of incorporating 
the stiffness parameter. The simplest variant is to multiply 
the corrections Ap by k £ [0...1]. However, for multiple 
iteration loops of the solver, the effect of k is non-linear. 
The remaining error for a single distance constraint after 
n s solver iterations is Ap(l — k)' h . To get a linear relation- 
ship we multiply the corrections not by k directly but by 
k' = 1 — (1 — fc)V"». With this transformation the error be- 
comes Ap(l — k')" s = Ap(l — k) and, thus, becomes linearly 
dependent on k and independent of n s as desired. However, 
the resulting material stiffness is still dependent on the time 
step of the simulation. Real time environments typically use 
fixed time steps in which case this dependency is not prob- 
lematic. 

3.4. Collision Detection and Response 

One advantage of the position based approach is how simply 
collision response can be realized. In line (8) of the simula- 
tion algorithm the M co u collision constraints are generated. 
While the first M constraints given by the object representa- 
tion are fixed throughout the simulation, the additional M co u 
constraints are generated from scratch at each time step. The 
number of collision constraints M co u varies and depends on 
the number of colliding vertices. Both, continuous and static 
collisions can be handled. For continuous collision handling, 
we test for each vertex i the ray x,- — > p,. If this ray enters an 
object, we compute the entry point q c and the surface normal 
n c at this position. An inequality constraint with constraint 
function C(p) = (p — q c ) • n c and stiffness k = 1 is added 
to the list of constraints. If the ray x,- — * p,- lies completely 
inside an object, continuous collision detection has failed at 
some point. In this case we fall back to static collision han- 
dling. We compute the surface point q. s which is closest to 
p; and the surface normal n s at this position. An inequality 
constraint with constraint function C(p) = (p — q. s ) • n s and 
stiffness k = 1 is added to the list of constraints. Collision 
constraint generation is done outside of the solver loop. This 
makes the simulation much faster. There are certain scenar- 
ios, however, where collisions can be missed if the solver 
works with a fixed collision constraint set. Fortunately, ac- 
cording to our experience, the artifacts are negligible. 

Friction and restitution can be handled by manipulating 
the velocities of colliding vertices in step (16) of the algo- 
rithm. The velocity of each vertex for which a collision con- 
straint has been generated is dampened perpendicular to the 
collision normal and reflected in the direction of the collision 
normal. 

The collision handling discussed above is only correct for 
collisions with static objects because no impulse is trans- 
ferred to the collision partners. Correct response for two dy- 



namic colliding objects can be achieved by simulating both 
objects with our simulator, i.e. the N vertices and M con- 
straints which are the input to our algorithm simply represent 
two or more independent objects. Then, if a point q of one 
objects moves through a triangle pi , P2, P3 of another object, 
we insert an inequality constraint with constraint function 

C(q,Pi,p 2 ,P3)=±(q-Pi)-[(P2-Pl) x (P3-Pi)] which 
keeps the point q on the correct side of the triangle. Since 
this constraint function is independent of rigid body modes, 
it will correctly conserve linear and angular momentum. 
Collision detection gets slightly more involved because the 
four vertices are represented by rays Xj — * p,-. Therefore the 
collision of a moving point against a moving triangle needs 
to be detected (see section about cloth self collision). 



3.5. Damping 

In line (6) of the simulation algorithm the velocities are 
dampened before they are used for the prediction of the 
new positions. Any form of damping can be used and many 
methods for damping have been proposed in the literature 
(see [NMK*05]). Here we propose a new method with some 
interesting properties: 

(1) Xcm = (Li x i m i)/(Li m i) 

(2) v cm = (£;Vim;)/CL-m;) 

(3) L = £ i r i x (m;V ( ) 

(4) I = 

(5) a = I- 1 L 

(6) forall vertices i 

(7) Av; = Vcm + tt> X r,- - V; 

(8) V/<-V/ + fcdampiiigAVi 

(9) endfor 

Here r, = x,- — x„„, f; is the 3 by 3 matrix with the property 
f,v = r,- x v, and damping £ [0. . . 1] is the damping coeffi- 
cient. In lines (l)-(5) we compute the global linear velocity 
x cm and angular velocity CO of the system. Lines (6)-(9) then 
only damp the individual deviations Av,- of the velocities v; 
from the global motion \ cm + fflxr,. Thus, in the extreme 
case damping = li om y the global motion survives and the 
set of vertices behaves like a rigid body. For arbitrary values 
of ^damping* the velocities are globally dampened but without 
influencing the global motion of the vertices. 



3.6. Attachments 

With the position based approach, attaching vertices to static 
or kinematic objects is quite simple. The position of the ver- 
tex is simply set to the static target position or updated at ev- 
ery time step to coincide with the position of the kinematic 
object. To make sure other constraints containing this vertex 
do not move it, its inverse mass w, is set to zero. 
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Figure 4: For bending resistance, the constraint function 
C(pi,P2,P3,P4) = arccos(ni • Ht) — <Po is used. The actual 
dihedral angle (p is measure as the angle between the nor- 
mals of the two triangles. 




Figure 5: Constraint function C(q,pi ,P2,P3) = (q — Pi) ■ 
n — h makes sure that q stays above the triangle Pi ,P2,P3 
by the the cloth thickness h. 



4. Cloth Simulation 

We have used the point based dynamics framework to im- 
plement a real time cloth simulator for games. In this section 
we will discuss cloth specific issues thereby giving concrete 
examples of the general concepts introduced in the previous 
section. 



4.1. Representation of Cloth 

Our cloth simulator accepts as input arbitrary triangle 
meshes. The only restriction we impose on the input mesh 
is that it represents a manifold, i.e. each edge is shared by at 
most two triangles. Each node of the mesh becomes a simu- 
lated vertex. The user provides a density p given in mass per 
area [kg/m 2 ]. The mass of a vertex is set to the sum of one 
third of the mass of each adjacent triangle. For each edge, 
we generate a stretching constraint with constraint function 

C st retch (P 1 , P2 ) = I P 1 - P2 1 - '() • 

stiffness k strelc i t and type equality. The scalar tg is the initial 
length of the edge and k stretc f, is a global parameter provided 
by the user. It defines the stretching stiffness of the cloth. For 
each pair of adjacent triangles (P1.p3.p2) and (pi,P2,P4) 
we generate a bending constraint with constraint function 

Cw(Pl>P2,P3,P4) = 

acos ( ( p2 ~ p ') X (P3~Pi) . (P2-P1) x (P4~Pi) \ _ 

\l(P2-Pl)x(P3-Pl)l l(P2 — Pi) X (P4-Pl)l/ ' 
stiffness kj, enc { and type equality. The scalar <po is the ini- 
tial dihedral angle between the two triangles and k\, enc j is a 
global user parameter defining the bending stiffness of the 
cloth (see Figure 4). The advantage of this bending term 
over adding a distance constraint between points P3 and P4 
or over the bending term proposed by [GHDS03] is that it is 
independent of stretching. This is because the term is inde- 
pendent of edge lengths. This way, the user can specify cloth 
with low stretching stiffness but high bending resistance for 
instance (see Figure 3). 

Eqns. (10) and (11) define the projection for the stretch- 
ing constraints. In the appendix A we derive the formulas to 
project the bending constraints. 



4.2. Collision with Rigid Bodies 

For collision handling with rigid bodies we proceed as de- 
scribed in Section 3.4. To get two-way interactions, we apply 
an impulse m/Ap,-/Af to the rigid body at the contact point, 
each time vertex i is projected due to collision with that 
body. Testing only cloth vertices for collisions is not enough 
because small rigid bodies can fall through large cloth tri- 
angles. Therefore, collisions of the convex corners of rigid 
bodies against the cloth triangles are also tested. 



4.3. Self Collision 

Assuming that the triangles all have about the same size, 
we use spatial hashing to find vertex triangle collisions 
[THM*03]. If a vertex q moves through a triangle pi, P2, 
P3, we use the constraint function 

(P2-Pl)x(P3-Pl) 



c(q,pi,P2,P3) = (q-pi)- TT^ Etvt^t r n "\i ~ h > 

l(P2-Pl) x (P3-Pl)l 

(12) 

where h is the cloth thickness (see Figure 5). If the vertex 
enters from below with respect to the triangle normal, the 
constraint function has to be 



C(q,Pi,P2,P3) = (q-pi) 



(P3-Pi) x (P2-P1) 

I (P3 — Pi ) X (P2-Pl) 



-h 
(13) 

to keep the vertex on the original side. Projecting these con- 
straints conserves linear and angular momentum which is es- 
sential for cloth self collision since it is an internal process. 
Figure 6 shows a rest state of a piece of cloth with self col- 
lisions. Testing continuous collisions is insufficient if cloth 
gets into a tangled state, so methods like the ones proposed 
by [BWK03] have to be applied. 

4.4. Cloth Balloons 

For closed triangle meshes, overpressure inside the mesh can 
easily be modeled (see Figure 7). We add an equality con- 
straint concerning all N vertices of the mesh with constraint 
function 

/"mangle* \ 
C(pi,...,p JV ) = ( £ (p ( , xp ( ,).p,j j - ^pressure Vq 

(14) 
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Figure 6: This folded configuration demonstrates stable self 
collision and response. 



and stiffness k = 1 to the set of constraints. Here t [ , tk and tK 
are the three indices of the vertices belonging to triangle i. 
The sum computes the actual volume of the closed mesh. It 
is compared against the original volume Vq times the over- 
pressure factor ^pressure- This constraint function yields the 
gradients 

V P, C = E (P,j x P,;)+ £ (P,iXP,;)+ £ (P^xPfj) 

r4=i r4=' v4=' 

(15) 

These gradients have to be scaled by the scaling factor given 
in Eq. (7) and weighted by the masses according to Eq. (9) 
to get the final projection offsets Ap,. 

5. Results 

We have integrated our method into Rocket [Rat04], a game- 
like environment for physics simulation. Various experi- 




Figure 7: Simulation of overpressure inside a character. 



ments have been carried out to analyze the characteristics 
and the performance of the proposed method. All test sce- 
narios presented in this section have been performed on a 
PC Pentium 4, 3 GHz. 

Independent Bending and Stretching. Our bending term 
only depends on the dihedral angle of adjacent triangles, not 
on edge lengths, so bending and stretching resistances can 
be chosen independently. Figure 3 shows a cloth bag with 
various stretching stiffnesses, first with bending resistance 
enabled and then disabled. As the top row shows, bending 
does not influence stretching resistance. 

Attachments with Two Way Interaction. We can sim- 
ulate both, one way and two way coupled attachment con- 
straints. The cloth stripes in Figure 8 are attached via one 
way constraints to the static rigid bodies at the top. In addi- 
tion, two way interaction is enabled between the stripes and 
the bottom rigid bodies. This configuration results in realis- 
tically looking swing and twist motions of the stripes. The 
scene features 6 rigid bodies and 3 pieces of cloth which are 
simulated and rendered with more than 380 fps. 
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Figure 8: Cloth stripes are attached via one way interaction 
to static rigid bodies at the top and via two way constraints 
to rigid bodies at the bottom. 



Real Time Self Collision. The piece of cloth shown in 
Figure 6 is composed of 1364 vertices and 2562 triangles. 
The simulation runs at 30 fps on average including self col- 
lision detection, collision handling and rendering. The effect 
of friction is shown in Figure 9 where the same piece of cloth 
is tumbling in a rotating barrel. 

Tearing and stability. Figure 10 shows a piece of cloth 
consisting of 4264 vertices and 8262 triangles that is torn 
open by an attached cube and finally ripped apart by a 
thrown ball. This scene is simulated and rendered with 47 fps 
on average. Tearing is simulated by a simple process: When- 
ever the stretching of an edge exceeds a specified threshold 
value, we select one of the edge's adjacent vertices. We then 
put a split plane through that vertex perpendicular to the edge 
direction and split the vertex. All triangles above the split 
plane are assigned to the original vertex while all triangles 
below are assigned to the duplicate. Our method remains sta- 
ble even in extreme situations as shown in Figure 1, a scene 
inspired by [ITF04]. An inflated character model is squeezed 
through rotating gears resulting in multiple constraints, col- 
lisions and self collisions acting on single cloth vertices. 

Complex Simulation Scenarios. The presented method 
is especially suited for complex simulation environments 
(see Figure 12). Despite the extensive interaction with an- 
imated characters and geometrically complex game levels, 
simulation and rendering of multiple pieces of cloth can still 
be done at interactive speed. 

6. Conclusions 

We have presented a position based dynamics framework 
that can handle general constraints formulated via constraint 
functions. With the position based approach it is possible to 
manipulate objects directly during the simulation. This sig- 
nificantly simplifies the handling of collisions, attachment 
constraints and explicit integration and it makes direct and 
immediate control of the animated scene possible. 

We have implemented a robust cloth simulator on top of 
this framework which provides features like two way inter- 
action of cloth with rigid bodies, cloth self collision and re- 
sponse and attachments of pieces of cloth to dynamic rigid 
bodies. 



7. Future Work 

A topic we have not treated in this paper is rigid body simu- 
lation. However, the approach we presented could quite eas- 
ily be extended to handle rigid objects as well. Instead of 
computing a set of linear and angular impulses for the res- 
olution of collisions as regular rigid body solvers typically 
do, movements and rotations would be applied to the bod- 
ies at the contact points and the linear and angular velocities 
would have to be adjusted accordingly after the solver has 
completed. 
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Figure 9: Influenced by collision, self collision and friction, a piece of cloth tumbles in a rotating barrel. 




Figure 10: A piece of cloth is torn open by an attached cube and ripped apart by a thrown ball. 




Figure 11: Three inflated characters experience multiple collisions and self collisions. 




Figure 12: Extensive interaction between pieces of cloth and an animated game character (left), a geometrically complex game 
level (middle) and hundreds of simulated plant leaves (right). 
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Appendix A: 

Gradient of the Normalized Cross Product 

Constraint functions often contain normalized cross products. To de- 
rive the projection corrections, the gradient of the constraint func- 
tion is needed. Therefore it is useful to know the gradient of a nor- 
malized cross product with respect to both arguments. Given the 
normalized cross product n = , the derivative with respect to 

the first vector is 
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Shorter and for both arguments we have 
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where p is the matrix with the property px = p x x. 



Bending Constraint Projection 

The constraint function for bending is C = arccos(rf) — <po, where 
d = ni ■ 112 = n^n 2 . Without loss of generality we set pi =0 and get 

for the normals ni = p^S? an( j n2 = P2 ><P4 win, d 

P2 X P3 IP: 
— J — - we get the following gradients: 
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Using the gradients of normalized cross products, first compute 
p 2 xn 2 + (n 1 xp 2 )d 
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Then the final correction is 
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